Eco-evolutionary dynamics of multigames with mutations

Most environments favor defection over cooperation due to natural selection. Nonetheless, the emergence of cooperation is omnipresent in many biological, social, and economic systems, quite contrary to the well-celebrated Darwinian theory of evolution. Much research has been devoted to better understanding how and why cooperation persists among self-interested individuals despite their competition for limited resources. Here we go beyond a single social dilemma since individuals usually encounter various social challenges. In particular, we propose and study a mathematical model incorporating both the prisoner’s dilemma and the snowdrift game. We further extend this model by considering ecological signatures like mutation and selfless one-sided contribution of altruist free space. The nonlinear evolutionary dynamics that results from these upgrades offer a broader range of equilibrium outcomes, and it also often favors cooperation over defection. With the help of analytical and numerical calculations, our theoretical model sheds light on the mechanisms that maintain biodiversity, and it helps to explain the evolution of social order in human societies.


Introduction
Numerous interdisciplinary researchers have long sought a way to understand how the tremendous biodiversity among species persists in nature despite the significant differences between them in terms of competitive capability. Darwinian evolution [1] always challenges the emergence of spontaneous social cooperation, as cooperators act entail an inherent cost for displaying altruism [2]. Self-interested defectors generally exploit these cooperators and hinder the maintenance of cooperation. This trivially raises the question of why cooperators act selflessly if only the fittest succeeds and how numerous forms of cooperative behavior sustain in nature ranging from a low-level microbial biological system to a high-level complex social system. We resort to the evolutionary game dynamical interaction given by a 2 × 2 pay-off matrix to solve this riddle. These classical games are capable of analyzing the economic and strategic decisions of rational individuals [3]. An unprecedented spectrum of researchers have already focused overwhelmingly on the evolutionary game theory [4][5][6] to seek out the mechanisms that sustain and promote cooperation [7][8][9][10][11][12][13][14].
We construct a mathematical model on evolutionary multigames [15][16][17][18][19][20][21][22][23][24][25][26][27] by adopting two simple two-player games, viz. the prisoner's dilemma [28] (PD) and snowdrift game [29,30] (SD), as metaphors for cooperation between unrelated rational individuals. Our choice of combining two distinct games is motivated by the fact that people may face different social complexities rather than a single one. People take care of problems based on their own perceptions. This real-life behavioral heterogeneity among individuals inspires us to examine the effect of multigames for understanding social diversity. Our theoretical approach is another effective way to explore the impact of the interaction between different cultural backgrounds on large-scale social ordering [31].
While defection is the ultimate profitable strategy in the classical PD game, irrespective of the co-player's decision, a player's best strategy in the SD game depends on the co-player's decision. Thus, we select two different games with two contrasting outcomes. The stable coexistence of both cooperators and defectors is the expected consequence in the SD game, resulting in the persistence of cooperative behavior. In comparison, cheaters are always encouraged to exploit the cooperative individuals in the set-up of the PD game. In our constructed model, unrelated individuals have their probabilistic selection to choose which game they want to play. We also consider the mutation [32][33][34][35][36][37][38][39][40][41][42] as an evolutionary mechanism that enables people to switch their respective strategies. We assume that it is hard to behave rationally under every circumstance. Hence, in our mathematical model, we incorporate a simplistic assumption that cooperators' and defectors' subpopulations can interchange their strategies with a fixed bidirectional mutation rate. To the best of our knowledge, the majority of the studies on multigames have been examined in the absence of mutation.
Apart from that, to reveal the interplay between ecological and evolutionary dynamics, we inspect the influence of free space in our mathematical model. In the context of complex systems, free space [43][44][45][46][47][48][49] proves to be a promising factor for offering significant consequences on diverse collective dynamics. We assume that each individual's birth rate depends not only on the respective average payoffs but are also proportional to the available free space. This free space will give them reproductive opportunities. The population of unrelated players gains a payoff due to the interaction of the multigames. We treat these payoffs as reproductive successes. Individuals with a higher payoff leave more offspring if sufficient empty spaces are available and are able to outcompete less successful ones. Free space will also contribute to the individual's average payoffs. We consider free space as an ecological variable donating selflessly to all players' fitness. However, free space never anticipates any benefit from others. This generous nature is surprising in a society of self-interested individuals. Still, free space provides reproductive opportunities to each player and loses their own identity to improve the fitness of other individuals. We consider such a selfless act of free space because helping others is a common practice among human beings and other animals. This tendency of unselfish concern for other people is a unique recipe for promoting cooperation and favors the survival of various species. The inclusion of such ecological dynamics into evolutionary multigames gives rise to a fascinating way to reveal the influence of eco-evolutionary dynamics [50][51][52][53][54][55][56][57][58][59][60] on decisionmaking and the evolution of cooperation.
The concurrence of ecological changes and rational players' evolution, depicted in the present article, prevails through the same time scale, rather than the typical assumption that ecological processes are much faster than evolutionary processes [61]. Our theoretical model uncovers how the combination of ecological dynamics and game dynamics is beneficial for maintaining cooperative behaviors under the influence of mutation and leads to the stable coexistence of interacting competitors. Earlier, Nag Chowdhury et al. [8,16] examined the influence of eco-evolutionary dynamics on cooperative behavior's emergence and evolution under a cooperation-supporting mechanism, viz. punishment. Nevertheless, our model exhibits the coexistence of cooperators and defectors in the absence of any such supporting mechanism. Moreover, individuals' birth rates in those Refs. [8,16] depend only on their respective average payoffs. Apart from the addition of free space's altruistic behavior in the payoff matrix, we furthermore include the importance of available free space on the birth rates by assuming their birth rates to be proportional to the available space. The rest of the article is structured as follows: In Sec. (2), we discuss how we construct our mathematical model in detail. The following section (3) deals with the existence, uniqueness, and boundedness of the model, along with extensive numerical simulations and discussions. The section (4) presents the summary of our findings. Following that, we provide a brief discussion in the last section (5) on the challenges in this work's context that need to be addressed and are worth studying in the future.

The model
To start with, we consider a simplistic assumption that each individual has two distinct choices, viz. (i) cooperation (C) and (ii) defection (D). Even they can play any of the two possible games (a) PD game and (b) SD game. They can adopt the PD game with probability p, and alternatively, they interact with other individuals by playing the SD game with the complimentary probability (1 − p). Both of these two-person games can be given by the following two payoff matrices A and B respectively, where in which the entries portray the payoff accumulated by the players in the left.
Here, R PD and R SD contemplate the reward for mutual cooperation among two players in the respective PD and SD games. Similarly, both unrelated players receive the punishment P PD and P SD for mutual defection in the games PD and SD, respectively. An exploited cooperator gains the sucker's payoff S PD and S SD , respectively in the PD and SD games when confronted by a defector. The mixed choice yields the defector temptations T PD and T SD to exploit a cooperator in the PD and SD games, respectively. The payoff ranking of these four-game parameters determines the two-person games. The conventional relative ordering for the PD game is T PD > R PD > P PD > S PD [8,16,31] and 2R PD > S PD + T PD [62]. Without loss of any generality, we choose T PD = β > 1, R PD = 1, P PD = η 2 (0, 1), and S PD = 0. Similarly, we choose T SD = β > 1, R SD = 1, S SD = 0, and P SD = −η 2 (−1, 0), maintaining the relative ordering T SD > R SD > S SD > P SD for the standard SD game [8,16,31]. Thus, the payoff matrices A and B become Note that, in both of these games, mutual cooperation leads to the payoff R PD and R SD , which is relatively higher than P PD and P SD , which one defector receives when playing with a defector. Thus, cooperation always promises higher income than defection if both the rational players choose the same strategy. The difference between these two games' relative ordering leads to a contrasting scenario. In the SD game, the interaction between the cooperator and defector always promises a better income in terms of payoff than the interaction between two defectors. A reverse reflection is observed in the case of the PD game thanks to the choice of such relative ranking of game parameters in the PD game. The interaction between two defectors in the PD game allows them to earn more than a cooperator encountering a defector. The switching between P and S in the relative ordering of both games thus produces noticeable unexpected consequences on the evolution of cooperation.
Since a player can decide which game they want to play, thus the final payoff matrix for the multigames looks like Here, p 2 [0, 1] is the probability of playing the PD game. Moreover, we consider free space F as an ecological variable, which contributes altruistically by helping others. Nevertheless, free space does not get any benefits by giving them reproductive opportunities. We incorporate this charitable role of free space by extending the 2 × 2 payoff matrix E to the 3 × 3 payoff matrix G as follows The matrix G clearly reveals that the free space never earns any payoff for their selfless charitable act; however, it contributes a positive payoff σ 1 and σ 2 to the cooperators and the defectors, respectively. Since most of the game parameters (not all) lies within the closed interval [0, 1], we assume, for the sake of feasible comparison, σ 1 and σ 2 both lies within the interval [0, 1]. When σ 1 and σ 2 are equal to zero, free space will not contribute anything to anyone. However, whenever σ 1 and σ 2 attain positive values, individuals gain an additional payoff from free space.
Inspired by the Malthusianism, we consider the following set of differential equations governing the changes in frequencies of cooperators and defectors as a function of time t Here, x and y are the normalized densities of cooperators and defectors, respectively. Let z be the available free space. Thus, we have Relation (2) assures that by studying x and y alone, one can easily capture the dynamics of the two strategies. We assume the birth rates of each individual depends crucially on the available free space as well as on their respective average fitness. Thus, we consider where f C and f D are average fitness of the cooperators and the defectors, respectively. The average payoff of cooperators and defectors can be determined using the payoff matrix G of the multigames, and the relation (2) as follows, Note that the average fitness of free space is This is expected as free space does not gain anything for its benevolent nature. For simplicity, we further assume that all individuals die at a uniform and constant mortality rate ξ 2 (0, 1]. Hence using the relations (3) and (4), our constructed model (1) transforms into Now, we introduce a constant probability μ as a rate with which each individual mutates from one strategy to the others, in a continuous manner, Relation (7) reflects that the mutation probability μ from the cooperators to the defectors is identical to the mutation rate from the defectors to the cooperators. So using the system (6), the well-mixed population under the influence of bidirectional mutation gives rise to the differential equations _ x ¼ x½ð1 À x À yÞfð1 À s 1 Þx À s 1 y þ s 1 g À x� þ mðy À xÞ; _ y ¼ y½ð1 À x À yÞfðb À s 2 Þx þ ðð2p À 1ÞZ À s 2 Þy þ s 2 g À x� þ mðx À yÞ: The system (8) contains seven different parameters. We summarize the necessary information about these parameters in (Table 1).

Existence, uniqueness and positive invariance
Before investigating the model (8) using numerical simulations, we first prove the positive invariance of the proposed system (8). Clearly, the functions on the right-hand side of the differential Eq (8) are continuously differentiable and at the same time, locally Lipschitz in the first quadrant of R � R. This ensures the existence and uniqueness of solutions for the model (8) with suitable initial values. Note that the initial conditions (x 0 , y 0 ) must lie within the domain [0, 1] × [0, 1] maintaining the inequality 0 � x 0 + y 0 � 1, as they represent the frequencies of the two strategies. To determine the positivity of the proposed model (8), we write the system as follows, where y are two integrable functions in the Riemannian sense. Solving Eq (9), we get where c 1 and c 2 are the integrating constants depending on the initial densities x 0 and y 0 . This proves that both x and y are non-negative. Now, to calculate the upper bound of x + y, we proceed as follows. The dynamical Eq (8) yield Since as per our previous analysis, x + y � 0 and the mortality rate ξ is always positive, we have ξ(x + y)�0. Thus, (11) reduces to Integrating both sides, we get where c 3 is the initial density dependent constant. Hence, we find that i.e., the overall species density x + y eventually remains bounded within the region [0, 1]. This boundedness within the closed interval [0, 1] allows us to relate the possible emerging dynamics of the system (8) to physically implementable scenarios with biological relevance. When the sum of the population density (x + y) is precisely one, the available free space is zero, as z = 1 − (x + y). i.e., there is no reproductive opportunity accessible to any individual in that situation with z = 0. When (x + y) = 0, then individually x = 0 and y = 0. Hence, all individuals die, and z

Coexistence of different stationary points depending on the initial conditions
Next, we point out the multistable dynamics of our model (8), resulting in the system's vulnerability to small perturbations. Initially, we set the parameter values at ξ = 0.38, β = 1.1, η = 0.85, p = 0.30, σ 1 = 0.36, σ 2 = 0.25 and μ = 0 in Fig 1(a). We vary the initial conditions (x 0 , y 0 ) within the interval [0, 1] × [0, 1] maintaining the inequality (x 0 + y 0 )2[0, 1]. We find that the dynamics switching between two stationary points, viz. E 0 = (0, 0) and E 1 = (x � , 0). We analytically calculate the stationary points of the system (8) in the absence of the mutation, i.e., with μ = 0. We trace out four different stationary points, 1. E 0 = (0, 0) reveals all individuals die. This point is locally stable when s 1 < x; and s 2 < x: ( 2. E 1 = (x � , 0) exhibits a society free from any defectors. Here, ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi The stability criteria is given by 3. E 2 = (0, y � ) represents that we have only left with defectors. Here, y � ¼ ð2pZÀ ZÀ 2s 2 Þ 2ð2pZÀ ZÀ s 2 Þ � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi indicates the coexistence equilibrium offering the survival of cooperation and defection simultaneously. Here, y �� ¼ ð1À bÀ s 1 þs 2 Þx �� þs 1 À s 2 2pZÀ Zþs 1 À s 2 , and x �� satisfies x �� and y �� both should lie within the interval (0, 1). The local stability yields the conditions for the stability of (x � , y � ) are Clearly, the chosen parameters satisfy the local stability criteria of both stationary points E 0 and E 1 . Note that E 1 leads to two different values for the selected values of the parameters, out of which (0.34751, 0) always remains locally stable, and (0.0899, 0) is unstable. Thus, we see the appearance of two different stationary points E 0 (shown by red points) and E 1 (displayed by violet points) in the basin of attraction portrayed in Fig 1(a). Such toggling between alternate stable stationary points is one of the generic features in some biological systems involving the fundamental processes of life [63][64][65] and in a few nonlinear dynamical systems [66][67][68][69]. The initial condition (x 0 , y 0 ) = (0, 0) always helps the system (8) to converge to the stationary point E 0 . The system is never able to give rise to the survival of any cooperators and defectors without the presence of any individual at the beginning. Since the chosen initial point (0, 0) itself is the stationary state, the system always stabilizes in the (0, 0) stationary point irrespective of the parameter values. Fig 1(b) is drawn with ξ = 0.15, β = 1.06, η = 0.85, p = 0.69, σ 1 = 0.75, σ 2 = 0.25, and μ = 0. Similarly, we find the system (8) in the absence of mutation (i.e., μ = 0) converges to E 0 for the single initial condition (x 0 , y 0 ) = (0, 0). It is anticipated that the choice of x 0 = 0 always leads to the cooperator-free steady state. The initial absence of cooperators in the mutation-free model will not entertain any exposure for the cooperators in the long run. The line of initial conditions x 0 = 0 and y 0 6 ¼ 0 produces the stationary state E 2 (sea blue points in Fig 1(b)). Apart from these two steady states, the mutation-free system also switches between E 1 (violet points) and E 3 (yellow points) depending on the suitable choice of initial conditions (See Fig 1(b)). Thus, for the same choices of parameters' values, the system flips between four alternate steady states depending on the initial densities of cooperators and defectors. Similarly, we find the system (8) with the parameters' values ξ = 0.15, β = 1.1, η = 0.85, p = 0.10, σ 1 = 0.10, σ 2 = 0.25 and μ = 0 converges to all these four stationary points. All population extincts for the initial conditions ranging from (0, 0) to (0.06, 0) (See red points in Fig 1(c)). The system (8) can be solved analytically with σ 1 = σ 2 = μ = 0 in the absence of defectors (i.e., y = 0) as follows, ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 4x À 1 p Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 4x À 1 p � : where c 4 is the integrating constant. Similarly, the system (8) can be solved analytically with σ 1 = σ 2 = μ = 0 in the absence of cooperators (i.e., x = 0) as follows 2Zð2p À 1Þ ; ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ð2p À 1ÞZ p tan À 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ð2p À 1ÞZ p ð2y À 1Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 4x À ð2p À 1ÞZ p ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where c 5 is the integrating constant. decide the asymptotic dynamics of the system. Other initial conditions with y 0 = 0 stabilize the dynamics in the defector-free state E 1 (violet points). We also track a fair portion of the basin in Fig 1(c), where the system (8) with μ = 0 converges to either E 2 (sea blue points) or E 3 (yellow points) depending on the choice of initial conditions.

> > > > > > > < > > > > > > > :
With increasing p, players are prone to play the PD game. Hence, we can not anticipate the survival of cooperators for large p. Thus, for large p, either the extinction equilibrium E 0 or the cooperator-free stationary point E 2 stabilizes in Fig 2(a), 2(c) and 2(e). However, the free-space induced benefits in Fig 2( With increasing β, the defectors are getting extra aid, and thus, we get the stabilization of the extinction equilibrium E 0 (red zone) in Fig 2(a). However, the presence of free space assistance converts that area into a coexistence zone (yellow zone) in Fig 2(b). The increment of p in both subfigures (a-b) permits only the existence of defectors, as people are playing mostly the PD game in such circumstances. Defectors are getting a favorable atmosphere in the PD game for our chosen PD game parameters' values. A similar sort of stabilization of cooperatorfree steady state (sea blue area) is observed in subfigures (c-d) for larger values of p. However, lower values of p provide the opportunity for playing the SD game. Hence, coexistence equilibrium (yellow region) stabilizes in Fig 2(c) and 2(d) for smaller p. Nevertheless, the inclusion of free space induced benefits helps to broaden the region of coexistence in the p − η parameter plane, as shown in Fig 2(d). The comparison between the Fig 2(e) and 2(f) suggests that appropriate non-zero values of σ 1 and σ 2 entertain the stabilization of a defector-free steady state (violet region). Hence, we can trace a fair portion of violet points in the p − ξ parameter space of the subfigure (f). However, too large a mortality rate reduces the opportunity of survivability of any individual, resulting in the stabilization of the extinction equilibrium E 0 . We focus on the effect of mortality rate ξ more elaborately in Fig 3. As expected, higher values of ξ constantly enlarge the chances of extinction. Thus, we observe a fair portion of the red region in Fig 3. Nevertheless, the amount of this red area is considerably lesser in the right column of Fig 3 compared to the left column. We introduce the non-zero values of σ 1 and σ 2 in the right column of this figure. These free space induced benefits encourage maintaining a defector-free society, as shown in Fig 3(b) and 3(d). We choose In the absence of free space induced benefits, we are unable to trace such defector-free regions in subfigures (a) and (c). The red region reflects the disappearance of all individuals for relatively high values of death rate ξ. Note that we keep the value of p fixed at 0.3. Thus, the system gets more opportunities to play the SD game, which generally encourages stabilizing the coexistence equilibrium. Therefore, we notice lower values of ξ will lead to the convergence towards the coexistence equilibrium E 3 (yellow points). β is kept fixed at 1.1 for the subfigures (a-b) and η is set at 0.85 for subfigures (c-d). We iterate the system (8) with μ = 0 for 20 × 10 5 iterations with fixed integration step size δt = 0.01 and fixed initial condition (x 0 , y 0 ) = (0.35, 0.35).
https://doi.org/10.1371/journal.pone.0272719.g003 σ 1 = 0.75 > 0.25 = σ 2 for Fig 3(b) and 3(d). i.e., the free space will provide an additional advantage for the cooperators, and thus depending on the other parameters, the defectors are vanished in the long run, as shown in Fig 3(b) and 3(d). In all of these subfigures of Fig 3, we trace a portion of yellow points depicting the survival of both cooperators as well as defectors simultaneously. We choose p = 0.3 in Fig 3. Thus, people get more chances to play the SD game, which facilitates the concurrence of both cooperation and defection. Thus, smaller values of ξ provide an opportunity to coexist for all strategies, which is observed in Fig 3 with  In Fig 4, we fix the value of σ 2 at 0.25 and examine the role of σ 1 . Since σ 1 represents the free space induced benefits towards the cooperators, thus enhancement of its (σ 1 ) value will help in . The color code represents the following: (i) red represents the extinct state, (ii) yellow portrays the co-existence state, (iii) violet displays the defector-free state, and (iv) sea blue depicts cooperator-free state, respectively. We have run the numerical simulations for 20 × 10 5 iterations for each point and store the final value for determining the final asymptotic state. Increment of σ 1 contributes more to the cooperators' payoff, and hence, we observe a defector-free region for higher values of σ 1 depending on the other parameters. We draw the stability curves for E 0 (the solid white line), E 1 (the blue dashed line), E 2 (the brown dotted line), and E 3 (the solid blue line). promoting cooperation in the society. Fig 4 reflects the same scenario. A defector-free society (violet region) is noticed in all these subfigures. The increment of temptation parameter β's value challenges the prevalence of cooperation. Thus, for a small σ 1 , we find a defector dominated society (the sea blue region in Fig 4(a)). σ 1 proves to be a cooperator facilitating parameter as we detect a wide range of yellow regions in Fig 4(a), where both cooperators and defectors can coexist. Similarly, irrespective of the choice of η in Fig 4(b), higher values of σ 1 provide all cooperators an extra benefit for survival, and thus, the stationary point E 1 (violet points) stabilizes. For the intermediate choice of σ 1 , the stationary point E 3 (yellow) yields the stable coexistence of both strategies. However, cooperators strive to keep in existence for the lower values of σ 1 , and we find the sea blue region of cooperator-free steady state. The parameter p indicates the probability of playing the PD game. Thus, p ! 1− always gives defectors a more favorable environment to survive. That's why we spot a sea blue portion in Fig 4(c) for higher values of p and σ 1 . When p is small, people are prone to play the SD game; and thus, we get the coexistence of both strategies (the yellow region in Fig 4(c)) for smaller p and σ 1 . Nevertheless, a larger value of σ 1 with a moderate value of p always provides a reasonable scope for . The color codes are as follows: (i) red represents the extinct state E 0 , (ii) yellow represents the co-existence state E 3 , (iii) violet represents the defector-free state E 1 , and (iv) sea blue represents cooperator free state E 2 , respectively. We plot the stability curves corresponding to each stationary point in each subfigure.
https://doi.org/10.1371/journal.pone.0272719.g005 the cooperators to survive, and we discover a healthy portion of the stationary point E 1 (the violet region) in Fig 4(c). A higher mortality rate ξ always results in the extinction of both cooperators and defectors, and we locate a huge red region in the σ 1 − ξ parameter plane of Fig  4(d). We find a very tiny sea blue region in Fig 4(d), where defectors are only able to survive. But, as expected, a higher value of σ 1 always promotes the cooperation strategy, and we obtain a violet zone of defector-free stationary point and a yellow region of interior equilibrium E 3 in Fig 4(d). Fig 5(a) shows that except for a smaller portion of the defector-free region (violet zone), the whole σ 2 − β parameter space produces the coexistence of both cooperators and defectors depending on other parameters' values. Despite the increment of σ 2 and β, the cooperators are able to survive along with the defectors due to our choice of other parameters' values. The σ 2 − η parameter space portrays that smaller values of σ 2 can not provide any benefit to the defectors, and stabilize the defector-free stationary point E 1 (violet region) irrespective of η's value. However, if σ 2 increases, it will yield a window of opportunity for the defectors to thrive. We observe a wide coexistence region (yellow region) and a small area of cooperator-free stationary point (sea blue region)in Fig 5(b). σ 2 indicates the free-space induced benefits towards the defectors. So, it is expected that larger values of σ 2 always enhance the chances of defectors' survivability. Thus, we notice a cooperator-free region (sea blue zone) in both Fig 5(c) and 5(d). Nevertheless, larger values of p enhance the chances of playing the PD game, where defectors get a favorable environment to survive. Thus, the cooperator-free sea blue region is found in Fig 5(c). We are also able to detect a small violet region of a defector-free environment in the σ 2 − p parameter space. However, we trace a healthy portion of coexistence (yellow) too in Fig 5(c) due to our chosen parameters' values.

The influence of bidirectional mutation
Now we investigate the influence of bidirectional mutation on the long-term behavior of the nonlinear differential Eq (8). Since every species can mutate into the other at a specific uniform rate μ 2 (0, 1], thus cooperators and defectors can not remain alive alone. Either all populations will die; otherwise, the dynamics will lead to the coexistence of all two species. Fig 4  already portrays the influence of free space-induced benefits on the cooperators in the absence of mutation. We scrutinize the impact of σ 1 under the influence of 50% mutation (i.e., μ = 0.5) in Fig 6. As σ 1 increases, the cooperators are getting a better environment for survival. We find a portion of stable coexistence equilibrium in each two-dimensional parameter space for large σ 1 in Fig 6. Thanks to the mutation, the cooperators can not live alone. The white line in Fig 6 is the stability curve corresponding to the extinction equilibrium. The local stability analysis fits almost exactly with the numerical simulations in Fig 6 with fixed initial condition (x 0 , y 0 ) = (0.35, 0.35). There are a few places where the stability analysis fails to predict the stabilization of the extinction stationary point (0, 0). This is mainly due to the multistable behavior of the proposed model 8.
As discussed, the increasing values of the parameters β, η, and p always provide the defectors a favorable environment to dominate the cooperators. However, beyond a critical value of σ 1 , both the cooperators and defectors can coexist, as shown in Fig 6(a)-6(c). The larger values of ξ always hinder the evolution of cooperators as well as of defectors. However, an intermediate choice of σ 1 − ξ favors the successful evolution of both strategies, as portrayed through Fig  6(d). The same feature is also noticeable in Fig 7. The complex evolutionary dynamics switch between two stationary points depending on the choices of parameters' values in Fig 7. For smaller values of σ 2 , the defectors are not getting enough advantages to survive, and thus the extinction equilibrium (0, 0) stabilizes in the red region of all subfigures of Fig 7. The choice of other parameters' values is also crucial for obtaining these stationary states. The parameter σ 2 benefits the defectors; hence the defectors can survive beyond a certain threshold of σ 2 . The employed 50% mutation rate helps to flow a certain fraction of defectors into cooperators, and we have a moderate portion of coexistence state (yellow region) in Fig 7. The observed results may vary for different choices of initial conditions, as the system is multistable. We plot the boundary separating solid white lines in all subfigures by analyzing the local stability analysis of the extinction equilibrium. Clearly, this stability curve agrees well with our numerical simulations, and the places, where they don't agree with the numerical simulations, is due to the multistable behavior of our proposed model 8. All the simulations are done by iterating for 20 × 10 5 times with fixed integrating step length δt = 0.01. The last point is gathered to finalize the asymptotic state. All the codes to generate these figures are publicly available at Ref. [70]. Fig 8 demonstrates the importance of the parameters σ 1 , σ 2 and ξ for the enhancement of cooperation under the presence of mutation. The larger values of σ 1 and σ 2 facilitate the evolution of at least one species, and the positive mutation rate μ 2 (0, 1] assures that species should mutate into the other. This mechanism will lead to the species' coexistence in a major portion (yellow region) of the two-dimensional parameter spaces represented in Fig 8(a) and 8(b). The results are further validated by plotting the stability curve (white solid lines) below which the  (Fig 4), which is drawn in the absence of mutation. Larger values of σ 1 always facilitate the maintenance of cooperation, and the bidirectional mutation reinforces the system's inherent tendency to flow from cooperators to defectors. The solid white line in each figure is the analytically computed stability curve corresponding to the extinction equilibrium.
https://doi.org/10.1371/journal.pone.0272719.g006 stationary point (0, 0) is locally stable. A higher mortality rate never entertains the evolution of both strategies; thus, we obtain the red region in Fig 8(c). This red region indicates the extinction equilibrium (0, 0). Once again, we plot the stability curve (solid white line) in Fig 8(c), above which both the species should be extinct as per our local stability analysis. All the subfigures are drawn with fixed initial condition (0.35, 0.35).
We inspect the influence of different parameters on the constructed model (8) in Fig 9. We keep fixed all parameters' values at ξ = 0.25, β = 1.1, μ = 0.5, p = 0.65, η = 0.85, σ 1 = 0.75 and σ 2 = 0.25, unless they are varied. Since free space provides additional advantage to the cooperators compared to the defectors as we choose σ 1 = 0.75 > σ 2 = 0.25, we have a defector-free society at least initially with μ = 0 in Fig 9(a). However, whenever the mutation rate μ becomes positive, each strategy can mutate into other. In this way, the density of the cooperators (red line) decreases, and the defectors' density (blue line) increases. Eventually, both densities almost become identical for μ ! 1 −. In of Fig 9(b), the vital role of p is investigated in the one-dimensional bifurcation diagram. As p ! 1 −, the defectors are getting the upper hand over the cooperators as p indicates the probability of playing the PD game. PD game always provides additional assistance to the defectors. All these results obtained in Fig 9 are consistent with the social dilemmas considered for constructing the model (8) . Fig 9(b) points out that lower values of p are always better for the maintenance of cooperation as small values of p indicate more rational people are playing the SD game and the SD game always favors the coexistence of both strategies. Increasing the temptation parameter β always uplifts the defectors' fraction y (blue line). Thus, one needs to choose the value of β wisely so that we obtain a moderate range of β in Fig 9(c) allowing the coexistence of both strategies. When β is small, we found the cooperators' density x (red line) dominates the defectors' fraction y (blue line). However for β > 1.45, y is larger than x. Note that the results may alter for a different choice of initial condition as the system 8 is multistable. We plot Fig 9 with fixed initial condition (0.35, 0.35) and the code to generate this figure is freely available at [70]. Increasing the parameter η 2 (0, 1) can provide extra benefits to the defectors. Thus, the rate of increment of y is slightly better Since the mutation rate μ is 50%, we can obtain only two stationary points of the model (8). Thus, initially, for a smaller choice of σ 1 , extinction of both species prevails in Fig 9(e). Nevertheless, with an increment of σ 1 , cooperators are getting further assistance from free space, and hence x becomes positive beyond σ 1 = 0.25. Since μ = 0.5, 50% of these cooperators mutate bidirectionally into the defectors, and thus, we have positive y too in that range of σ 1 . Since σ 2 = ξ = 0.25, thus we get the stabilization of (0, 0) within the interval σ 1 2 [0, 0.25]. In fact, the Effects of various parameters in spreading of cooperative behavior. Subfigure (a) reveals two contrasting scenarios, where the cooperators decrease with the increment of the mutation rate μ. However, the defectors' fraction gradually increases with increasing μ. In the absence of mutation (i.e., for μ = 0), all defectors are extinct as the parameters values chosen to produce this subfigure is ξ = 0.25, β = 1.1, η = 0.85, p = 0.65, σ 1 = 0.75 and σ 2 = 0.25. Thus, the free space induced benefit towards cooperator σ 1 is greater than that of defector σ 2 . Subfigure (b) shows that both species' densities gradually increase with the probability p of playing the PD game. As the probability of playing the PD game increases, the increment of the fraction of defectors (blue line) gets over in margin with respect to the fraction of the cooperators (red line) after p = 0.85 depending on other parameters' values. This fact is obvious from the view that a higher possibility of playing the PD game gets a more beneficial ambiance for defectors. Subfigure (c) indicates the same phenomena as the value of the temptation payoff parameter β increases, both the fractions escalate. Whereas at β = 1.45, both these fractions get the same value. After that value of β, the value of the fraction of defector (blue line) becomes more than that of the cooperators (red line). Subfigure (d) also highlights the increment of both the fractions as the punishment payoff parameter η increases. But up to the highest received payoff value of η, the value of the fraction of defectors can never exceed the cooperators' fraction because the free space provides more benefit to the cooperators than the defectors. Subfigure (e) demonstrates that when σ 1 2 [0, 0.25], both these fractions are extinct, as the mortality rate ξ has more value than these free space-induced benefits. But as the value of σ 1 increases, the cooperator gets more advantages than defectors. Subfigure (f) shows that both the fractions increase with the enhancement of σ 2 . Up to σ 2 = 0.5, the cooperators' density acquires more value than y, as the value of σ 1 = 0.75 is taken sufficiently high. After σ 2 = 0.5, y becomes higher than x. All the subfigures are drawn by iterating the system (8)  cooperators' density x is slightly better than that of y for the higher values of σ 1 in Fig 9(e). The choice of the initial condition and other parameters' values are vital in obtaining all these results. Similar observation can be found in Fig 9(f), where we examine the role of σ 2 in our proposed model (8). With increasing σ 2 , the defectors will be benefited from free space, and hence, y (blue line) dominates x (red line) for larger values of σ 2 . But since σ 1 is taken sufficiently large in this Fig 9(f), we initially have a small portion for smaller values of σ 2 where x > y. Fig 9 discloses suitable choices of all parameters' values not only entertain the coexistence of both strategies but also may promote the evolution of cooperation.

Concluding remarks
Our mutation-induced model provides widespread coexistence of both strategies under suitable choices of parameters' values. Such stable persistence of ecological communities strengthens the theory of concurrency. Our thorough analysis with several numerical simulations enhances our understanding of the mechanisms that drive the survival of cooperative behavior in the multigames consisting of both the Prisoner's Dilemma and the Snowdrift games. The inclusion of altruistic behavior of free space along with the mutation in our proposed evolutionary model seems to be a natural course of action as observed in many realistic settings. Our findings clearly demonstrate that there exists an optimal probability of playing each game so that people are more likely to cooperate in such circumstances. Throughout the study, we have pointed out the positive impact of diverse factors (parameters) on significant improvement of cooperation. We have shown that the selfless contribution of free space promotes the coexistence of all strategies efficiently, even in the absence of mutation.
In summary, our research indicates the proposed model (8) may possess four different stationary points in the absence of mutation. The exciting feature of this model is that the system never allows settling into a cooperator-free state in the lack of free space-induced benefits and mutation rate. This precise result is consistent with the chosen games as both the PD and SD games never encourage a cooperator-free society in the usual scenario. The merging of two games with different outcomes provides a more realistic representation of the concept of opinion formation. However, the numerical results presented here are highly sensitive to the variation of initial conditions. The positive invariance and boundedness of the model are analyzed too in this article. We have shown the viable choice of bidirectional mutation allows the system to switch between only two stationary states. Either all people will die, or both the strategies coexist in the eco-evolutionary model with mutation. This is an interesting angle of our research as our model brings forth stable biodiversity in the form of a heterogeneous population (mixed cooperator-defector state). Note that we have only focused on the equilibria of the corresponding dynamical system (8) throughout the article so that we can relate those stationary states from the game-theoretical point of view. Our simple model with mutation sheds light on how cooperation emerges in a complex society. The presented insightful results attest that altruistic behavior and mutation are advantageous for the spontaneous maintenance of biodiversity. We believe the presented results may help us understand the mechanism behind the coexistence of competing species through the co-evolution of both strategies.

Discussion: Limitation & future perspectives
Lastly, we discuss challenges related to our work that merits further investigation. Even though it is almost impossible to incorporate all complex, realistic relationships among social creatures using minimal modeling, it is nevertheless absolutely essential to investigate such models using elementary mathematical principles. In fact, the literature already provides a number of excellent models that offer a stimulating starting point for exploring many practical situations.
Chen et al. [71] proposed an elegant model aimed at understanding the co-evolutionary outcome of the strategies of the traffic management department, drivers, and pedestrians. Reference [72] proposes a hybrid machine learning model to predict traffic accidents. A novel autonomous system without any stationary points may offer hidden attractors, as revealed in Ref. [73]. It is almost impossible to incorporate even a partial list of relevant references to emphasize the importance of minimal modeling here. Motivated by all these works, we present a new model that can offer several valuable insights into the evolution of cooperation and altruism. However, our model is far from perfect. Instead of considering two-person interactions, the introduction of group interactions may yield a deeper understanding of decisionmaking. Interdependent networks may be an excellent choice to study the evolution of cooperation, as reflected through Refs. [20,31,62,74]. This remains a promising future research generalization, including interdependent networks which may yield several complex dynamical behaviors other than stable equilibrium states.
One should observe that our approach works absolutely fine with other social dilemmas.
For any two two-person games of the form R S T P � � , the system constructed with the help of our policy always remains bounded within the closed interval [0, 1]. However, the stationary points alter due to the change in the mathematical model. We have verified this boundedness by considering the snowdrift and the harmony game [75,76] (results are not shown here). In fact, on the generality of our results, we can comment that a similar impact of mutation and free space can also be expected for other two-person games. It will be interesting to examine how the results may vary with other games like public goods games, rock paper scissor, etc. Besides introducing other realistic scenarios like apology and forgiveness [77], intention recognition [78], delay [79] etc., one may advance our understanding of the roots of cooperation in social and biological systems. Moreover, the impact of higher-order interactions [80] among the agents on the coevolution [81] of cooperation and synchronization [82][83][84] in a coupled network remains largely unexplored. The interdisciplinary researchers of complex systems can pay attention to this exciting topic for future research.